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Evaluation of a Distributed Catchment Scale 
Water Balance Model 

Peter A. Troch . 1 Marco Mancini,- Claudio PANIcoNL , • 4 and Eric F. Wood" 

The validity of some of the simplifying assumptions in a conceptual water balance model is 
investigated by comparing simulation results from the conceptual model with simulation results from 
a three-dimensional physically based numerical model and with field observations. We examine, in 
particular, assumptions and simplifications related to water table dynamics, vertical soil moisture and 
pressure head distributions, and subsurface flow contributions to stream discharge. The conceptual 
model relics on a topographic index to predict saturation excess runoff and on Philip s infiltration 
equation to predict infiltration excess runoff. The numerical model solves the three-dimensional 
Richards equation describing flow in variably saturated porous media, and handles seepage face 
boundaries, infiltration excess and saturation excess runoff production, and soil driven and atmo- 
sphere driven surface fluxes. The study catchments (a 7.2-km* catchment and a 0.64-km~ subcatch- 
ment) are located in the North Appalachian ridge and valley region of eastern Pennsylvania. 

Hydrologic data collected during the machydro 90 field experiment are used to calibrate the models 
and to evaluate simulation results. It is found that water table dynamics as predicted by the conceptual 
model are close to the observations in a shallow water well and therefore, that a linear relationship 
between a topographic index and the local water table depth is found to be a reasonable assumption 
for catchment scale modeling. However, the hydraulic equilibrium assumption is not valid for the 
upper 100 cm layer of the unsaturated zone and a conceptual model that incorporates a root zone is 
suggested. Furthermore, theoretical subsurface flow characteristics from the conceptual model are 
found to be different from field observations, numerical simulation results, and theoretical baseflow 
recession characteristics based on Boussinesq’s groundwater equation. 


1. Introduction 

Scale effects on hydrologic response have received much 
attention during the past decade. One of the central issues 
concerns the extent to which small scale details of the 
relevant physical processes should be represented in hydro- 
logic models applicable over a wide range of scales. Past 
research efforts have resulted in the formulation of scale- 
related concepts that summarize basin characteristics. Ex- 
amples can be found in the works by Gupta and Waymire 
[1983] and Mesa and Mifflin [1986], who used the network 
width function to describe channel flow routing; Gupta and 
Waymire [1989] and Tarboton et a!. [1989], who developed a 
stochastic theory to describe spatial variability in link 
heights in river networks; Seven and Kirkby [1979] and 
Sivapalan et al. [1987], who used the topographic index 
distribution of a drainage basin to model saturation excess 
runoff; and Wood et al. [1988] and Wood et al. [1990], who 
introduced the concept of the representative elementary area 
(REA) as the scale at which continuum assumptions con- 
cerning small scale heterogeneity hold. 

Based on these ideas, conceptual models have been de- 
veloped to describe and explain hydrologic response at the 
basin scale. What these models have in common is their 


'Laboratory of Hydrology and Water Management, Ghent Uni- 
versity, Gent, Belgium. 

2 Dipartimento di lngegneria Idraulica, Ambientale e del Rileva- 
mento, Politecnico di Milano, Milan, Italy. 

^Dipartimento di Metodi e Modelli Matematici, Universita di 
Padova, Padua. Italy. 

4 Now at CRS4, Cagliari, Italy. 

s Water Resources Program, Department of Civil Engineering and 
Operations Research. Princeton University, Princeton, New Jersey. 

Copyright 1993 by the American Geophysical Union. 

Paper number 93WR00398. 

0043- 1 397 93/93 WR-0039HS05. 00 


representation of physical process dynamics using analyti- 
cally tractable solutions to governing equations. To obtain 
such closed form solutions, it is in general necessary to 
consider simplifications in the mathematical formulation of 
the processes, such as linearizing the governing equations or 
assuming a restricted set of initial and boundary conditions. 
Moreover, different and often contrasting analytical solu- 
tions are used to represent the various hydrologic processes 
of interest. It is then hoped that the hydrologic model will be 
capable of simulating some average response of the basin. In 
this spirit, calibration and validation of conceptual models is 
often based on lumped hydrologic fluxes such as total 
discharge at the catchment outlet. However, current usage 
of conceptual models is not restricted to prediction of 
average hydrologic behavior. In environmental studies, for 
instance, accurate prediction of the temporal and spatial 
structure of each component of the hydrologic cycle has 
become increasingly important. 

It is the purpose of this paper to evaluate some of the 
simplifying assumptions in a distributed basin scale water 
balance model. The conceptual model under study is de- 
scribed in detail by Famiglietti et al. [1992]. The model relies 
on a topography-soil index to predict saturation excess 
runoff [Be ven and Kirkby . 1979; Sivapalan et al 1987] and 
on Philip's equation [ Brutsaert , 1976; Philip , 1957a, b] to 
predict infiltration excess runoff. In this paper, we examine 
conceptualizations concerning the spatial distribution and 
temporal evolution of the water table depth, the distribution 
of soil moisture and pressure head in the unsaturated zone, 
and the characteristics of base flow' recession. The topogra- 
phy-soil index, an index of hydrologic similarity \Wond et 
al. , 1990], is used to compute the local depth to the water 
table. In the model this water table depth will determine the 
storage capacity and hydraulic properties in the unsaturated 
zone, and it therefore plays a central role in the generation of 
surface runoff. The local storage deficit is calculated based 
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on a hydraulic equilibrium assumption for pressure head, 
from which the vertical distribution of soil moisture is 
obtained. This hydraulic equilibrium assumption also yields 
the surface soil moisture value which in turn is used to 
calculate the parameters in Philip's infiltration equation. It is 
further assumed that during a rainfall event the average 
water table depth remains constant. Immediately after ces- 
sation of rainfall the change in groundwater Storage is 
proportional to the infiltrated and drained volume, and this 
information is used to update the water table depth. During 
an interstorm period the water table depth is a function of 
evaporation and subsurface flow', and is updated at every 
time step. In the model subsurface contributions to total 
streamflow depend on an exponential relationship between 
base flow and groundwater storage. 

To investigate the assumptions described above we com- 
pare simulation output from different components of the 
conceptual model with results obtained from a three- 
dimensional numerical model and with field observations. 
The numerical model used in this study is described by 
Paniconi and Wood [1993] and is based on the three- 
dimensional transient Richards equation. This model han- 
dles evaporation and precipitation inputs at the catchment 
surface and predicts saturation excess and infiltration excess 
runoff generation. The evaluation of the conceptual model is 
carried out at two basin scales: for a 7.2-knrp catchment and 
for a 0.64-kirr subcatchment. These study catchments 
(WE-38 and WD-38 Mahantango Creek) are located in the 
North Appalachian ridge and valley region of eastern Penn- 
sylvania. Detailed hydrologic information for both catch- 
ments was collected during the 12-day Multisensor Airborne 
Campaign 1990 (machydro 90). One of the objectives of this 
campaign is to apply the multisensor data in detailed water 
balance studies of the regulating effect of soil moisture on 
hydrologic processes. 

Numerical models based on the Richards equation have 
been used in the past to simulate hillslope and catchment 
scale hydrologic processes [e.g.. Freeze , 1971; Smith and 
Hebhert , 1983; Binley et al , 1989]. These physically based 
models have proved useful in evaluating the underlying 
assumptions in conceptual models. Examples can be found 
in the works by Gan and Burges [1990a, b ] (rainfall-runoff 
models on small hypothetical catchments), Sloan and Moore 
[1984] (one- and two-dimensional subsurface storm flow 
models), and Ibrahim and Brutsaert [1968] and Reeves and 
Miller [1975] (one-dimensional infiltration models to test the 
time compression approximation). The numerical model 
used in this study is based on the following assumptions and 
limitations: flow is laminar and isothermal, inertial forces 
and chemical gradients are neglected, and the air phase is 
continuous and at atmospheric pressure. In addition, the 
model does not account for hysteresis and it is assumed that 
the porous medium is isotropic. Finally, we consider only 
flow within the soil matrix, neglecting flow- through macro- 
pores. 

In the following section we present the main features of 
the conceptual model and the numerical model. Both models 
can handle spatially and temporally variable inputs and are 
designed to take advantage of digital elevation data bases 
and of information extracted from these data bases by 
topographic analysis. Mechanisms of streamflow generation 
are discussed. Geophysical characteristics of the study 
catchments arc given in section 3.1. The Mahantango Creek 


watershed is characterized by long even-crested ridges, 
which alternate with broad, rolling valleys. The catchment 
experiences a humid climate w'iih approximately 1000 mm of 
annual precipitation, distributed uniformly throughout the 
year. The data and parameters used in the two models are 
described in sections 3.2 and 3.3, with additional discussion 
in section 3.4 of initial conditions and the estimation of an 
effective depth to the water table. The calibration of both 
models is discussed in section 3.5 and is based on a compar- 
ison betw-een observed and simulated runoff volume. 

In the conceptual model the local water table depth is 
assumed to be a linear function of the topography-soil index. 
This hypothesis is tested in section 4.1 by comparing the 
distribution of the topography-soil index with the distribu- 
tions of the water table depth generated by the numerical 
model for subcatchment WD-38 at the end of a 12-day 
simulation period. To minimize the effect of initial condi- 
tions, longer test runs with the numerical model are also 
performed. In section 4.2 assumptions in the conceptual 
model about the temporal evolution of the local and average 
water table depth are studied. Piezometric observations 
from a shallow and a deep well are compared with simulation 
results from both the conceptual and numerical models. In 
the following section the hydraulic equilibrium assumption 
for the vertical distribution of soil moisture is tested by 
generating moisture profiles at fixed points along a hillslope 
situated in subcatchment WD-38 using the numerical model. 
In section 4.4 the characteristics of base flow recession for 
the conceptual and numerical models are compared with 
observations and with analytical solutions to Boussinesq’s 
hydraulic groundwater equation. 

2. Description of the Models 

2. 1 . Conceptual Water Balance Model 

The conceptual water balance model used in this study is 
developed by Sivapalan et al. [1987] and Famiglietti et al. 
[1992]. We present here a brief description of the basic 
equations of the model. The water balance for a catchment 
with drainage area A is given by 


Ae = Ap 



dS g 

dt 


(1) 


where / is time, e is the catchment average evaporation rate, 
p is the lumped rainfall intensity, Q is the streamflow at the 
basin outlet, 5 W is the volume of water stored in the 
unsaturated zone, and S g is the groundwater storage. Other 
contributions to the water balance, such as deep percolation 
and the storage in surface water bodies, are assumed to be of 
minor importance and are neglected. 

2.1.1. Streamflow generation. The conceptual model 
predicts soil saturation and its relationship to both saturation 
excess and infiltration excess surface runoff generation using 
spatially variable soil and topographic data. The initial 
storage capacity in the unsaturated zone is a function of the 
depth to the w'ater table. The local depth to the water table 
Zi can be computed as {Sivapalan et al 1987) 

where 
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z areal average of e ( ; 

Tj local value of the transmissivity coefficient / = 

K s saturated hydraulic conductivity at the surface; 

/ describes the exponential decay ot hydraulic 
conductivity with depth; 

In (T e ) expected value of In ( 7',); 

A expected value of the topographic index In (u/tan 

0); 

a local draining area per unit contour length; 
p land surface slope angle. 

The derivation of (2) is based on the assumption that the 
water table is nearly parallel to the soil surface and that the 
saturated hydraulic conductivity K x declines exponentially 
with depth z 0 t according to 

K s iz) = K So exp i-fz d ) (3) 


Assuming hydraulic equilibrium for the vertical pressure 
distribution, the local storage capacity 5/ can be expressed 
in terms of z, and given soil parameters as 
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where 6 S is the saturation moisture content and 0 r is the 
residual moisture content. The air entry pressure head value 
i// c and the pore size distribution index B are the parameters 
in the soil moisture characteristic relationship formulated by 
Brooks and Corey [1964]: 

0(xP) = 6 r + (B s ~ * =£ 

\* } (5) 

e(iif) = e s ^ 


K{ * ) = Ks \i j 

K(\p) - K 5 ^ ^ i P c 


where 0 is volumetric moisture content, ip is pressure head, 
and K is hydraulic conductivity. Whenever the cumulative 
infiltration volume exceeds the local soil moisture storage 
capacity, saturation excess runoff is generated. 

The infiltration excess component of the model applies the 
time compression approximation to Philip’s equation to 
compute the local infiltration rate g { under local erratic 
rainfall /?, [Famiglietti et aL , 1992]: 

9i = min [$r*(G), p,] (7) 


The model assumes an exponential relationship between 
base flow and groundwater storage such that the base flow 
contribution Q h to total streamffow is given by 

G/> - Go exp (-/:) < l >> 

where Go represents subsurface flow when c = 0. During a 
rainstorm the average depth to the water table ; is kept fixed 
and redistribution of soil moisture is neglected. After cessa- 
tion of rainfall a new value tor z is calculated taking into 
account the total amount of water infiltrated and drained. 

The total streamffow for the catchment is obtained by 
summing the contributions to surface runoff from saturation 
excess, infiltration excess, and subsurface flow. In this 
version of the model overland flow and channel flow' routing 
are not considered. 

2.1.2. Evaporation . Saturated areas are allowed to 
evaporate at the potential rate c p (atmosphere controlled 
stage). As the soil near the surface dries out the moisture 
delivery rate is limited by the properties of the soil profile 
(soil controlled stage). The model uses an analytical solution 
to the one-dimensional desorption problem and again applies 
the time compression approximation which results in the 
following expression: 

e a = min [**(£„), e p ] (10) 

in which e a is the actual evaporation rate and >* represents 
the evaporation capacity as a function of cumulative evapo- 
ration E a : 

el(E a ) = S;/2E fl <*1> 

where S c is the desorptivity, which varies with soil moisture 
content. Transpiration is not considered in this version of the 
model. 

During interstorm periods the water table is updated at 
each time step, taking into account evaporative and drainage 
losses. 


2.2. Numerical Model 

The numerical catchment simulation model is presented in 
the works by Paniconi [1992] and Paniconi and Wood [1993]. 
We will highlight the main features of the model. The 
three-dimensional Richards equation w ith pressure head ip as 
the dependent variable can be written as 

d ip 

S(iP) — = V - (K s K r (iP)V(iP + Z)) (12) 


in which G is cumulative infiltration and the local infiltration 
capacity is calculated as 


g*(G\ - 1 + 


4rA\ G 


s; 


1/2 


- 1 


( 8 ) 


where S r is sorptivity and cA Si accounts for the effect of 
gravity on the infiltration rate. The parameter c ranges from 
0.5 to 1.0 and depends on the saturation of the soil profile. 
Analytical expressions for c and .S r are given in the work by 
Siva pal an and Wood [1986]. Infiltration excess runoff is 
generated on those parts of the catchment where p t > g*. 


where z is the vertical coordinate, positive upward, and the 
hydraulic conductivity K is expressed as a product of the 
conductivity at saturation and the relative conductivity A r , 
An extension of the van Genuchten characteristic equations 
[van Genuchten and Nielsen . 1985] is used to describe the 
nonlinear dependencies of ft, K r , and specific moisture 
capacity S on the pressure head [ Paniconi . 1992); 

0{ip) = 6 r + (0, - 0,)[1 ^ pV m iP^iPv 
6 UP) = ft r + (0, ~ ft r )[l + P 0 ]"'" + SsUP ” </'«) «/' - 
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T ABLE I. Soil Parameters lor WE -38 Mahantarij>o Creek 

Water Retention Parameters 


Soil 

Number 

Soil Name 

Soil 

Texture 

71 

Albrights 

silt loam 

57 

Alvira 

silt loam 

8 

Basher 

silt loam 

145 

Berks 

silt loam 

45 

Calvin 

silt loam 

73 

Conyngham 

silt loam 

32 

Dekalb 

sandy loam 

54 

Hartleton 

silt loam 

149 

Klinesvil 

silt loam 

75 

Laidig 

gravel loam 

66 

Leek Kill 

silt loam 

69 

Meckesville 

loam 

70 

Meckesville 

stony loam 

38 

Shelmadine 

silt loam 

47 

Weickert 

silt loam 


van 

Brooks-Corey Genuchten 

^ 

m 'hour 0, 0 T <//, . m B ,i, . m n 


0.036 

0.501 

0.015 

0.036 

0.501 

0.015 

0.036 

0.501 

0.015 

0.073 

0.501 

0.015 

0.057 

0.501 

0.015 

0.036 

0.501 

0.015 

0.090 

0.450 

0.041 

0.057 

0.501 

0.015 

0.090 

0.501 

0.015 

0.090 

0.460 

0.027 

0.057 

0.501 

0.015 

0.036 

0.501 

0.015 

0.090 

0.501 

0.015 

0.036 

0.501 

0.015 

0.073 

0.501 

0.015 


0.21 

0.211 

0.43 

1.29 

0.2! 

0.211 

0.43 

1.29 

0.21 

0.211 

0.43 

1.29 

0.21 

0.211 

0.43 

1.29 

0.21 

0.21! 

0.43 

1.29 

0.21 

0.211 

0.43 

1.29 

0.15 

0.322 

0.24 

1.38 

0.21 

0.211 

0.43 

1.29 

0.21 

0.211 

0.43 

1.29 

0.11 

0.220 

0.17 

1.25 

0.21 

0.211 

0.43 

1.29 

0.21 

0.211 

0.43 

1.29 

0.21 

0.21! 

0.43 

1.29 

0.21 

0.211 

0.43 

1.29 

0.21 

0.211 

0.43 

1.29 


5 ( 0 ) = —= (n ~ l)ie >~ 6 ^\ n ' 1 

d* \<l' s \ n (l + P ) m * 1 


^ <Ao 


dO 


(14) 


K r {W = {\ + /3) ~ 5m/2 [(i + p)'»-p*]2 ^ < 0 


A" r (0) =1 0 > 0 

where 5* is the specific storage, m = I - 1/n, /3 - (0/0J /l , 
/3o s is a fitting parameter, and n can be 

interpreted as a pore size distribution index; 0 O is a conti- 
nuity parameter which is calculated from (14) given S s . The 
exponential relationship (3) is used to model vertical heter- 
ogeneity of saturated hydraulic conductivity. Hysteresis 
effects on moisture redistribution in the soil profile are not 
taken into account. We remark here that different soil 
moisture characteristic relationships are used in the concep- 
tual model (equations (5) and (6)) and the numerical model 
(equations (13) and (15)). Any discrepancy this may cause is 
minimized by fitting both sets of equations to the same field 
observations (see Table 1). Since the soils in Mahantango 
Creek show a wide pore size distribution, the van Genuchten 
water retention equation fitted the observed data better than 
the Brooks-Corey equation. Track et al. f 1 993] and P. A. 
Troch et al. (Hydrologic controls of large floods in a small 
basin: North Appalachian case study, submitted to the 
Journal of Hydrology , 1993) have developed a more general 
formulation of the conceptual model based on the van 
Genuchten characteristic equations. The effect of the use of 
different soil characteristic curves in catchment modeling is 
under investigation. 

A finite element Galerkin discretization in space and a 
finite difference discretization of the time derivative term is 
used to solve (12) numerically. The resulting system of 
nonlinear equations is linearized using either Picard or 
Newton iteration [Paniconi et al. . 1991]. 

I he catchment simulation model comprises two programs; 


a grid generator w-hich constructs the finite element mesh 
and initializes various parameters, and the actual simulation 
program which numerically solves the three-dimensional 
Richards equation over a specified time period for a given set 
of boundary and initial conditions. 

2.2.1. Streamftow generation . The potential inflows to 
the model consist of precipitation (positive) and evaporation 
(negative) flux inputs at the catchment surface. The actual 
(simulated) inflows are determined according to the type of 
boundary condition imposed, and during a simulation the 
model automatically adjusts this boundary condition accord- 
ing to changes in pressure head and flux values at the 
surface. When the potential flux is positive, the difference 
between potential and actual soil inflow is the total runoff. 
Surface runoff is produced when the surface becomes satu- 
rated, either due to a rising water table (saturation excess 
mechanism) or to the infiltration capacity of the soil falling 
below the rainfall rate (infiltration excess mechanism). In 
both cases the boundary condition at the point on the surface 
where saturation occurs switches from a Neumann type 
(atmosphere controlled inflow') to a Dirichlet type (soil 
controlled inflow). Subsurface runoff in the model is pro- 
duced at seepage faces or when subsurface water exits the 
soil matrix from a saturated region on the surface (return 
flow). Overland flow (surface runoff and return flow) gener- 
ated at a point on the catchment is routed to the stream using 
a time delay determined from the overland flow velocity 
(assumed constant for the catchment) and the shortest 
distance from the point to the stream. In this version of the 
model channel flow is not considered. 

2.2.2. Evaporation . When evaporation is atmosphere 
controlled and the pressure head at a point on the surface 
becomes smaller than the “air dry” pressure head value 0 mjn 
\HilleL 1980, p. 121J, the boundary condition at that point is 
switched from specified flux (Neumann) to constant head 
(Dirichlet) and thus the evaporation process becomes soil 
controlled. The boundary condition switches back to a 
Neumann type when the magnitude of the computed flux 
across the soil surface exceeds the magnitude of the poten- 
tial evaporation rate, or when a rainfall event begins. Root 
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extraction is not included in the numerical model and there- 
fore the model cannot simulate evaporation from vegetated 
surfaces when the actual rate falls below the potential rate. 

3. Description of the Catchments 
and Model Inputs 

3.1. Mahantango Creek Experimental Catchment 

The Mahantango Creek watershed (approximately 420 
km 2 ) is characterized by long, even-crested ridges 300-400 
m in elevation, which alternate with broad, rolling valleys, 
150-300 m in elevation The mountain ridges run in the 
NE-SW direction. The ridges are underlain by erosion- 
resistant sandstones, orthoquartzite, and conglomerate. The 
valleys are underlain by less erosion-resistant shales, silt- 
stones, and sandstones [ Urban , 1977]. The catchments used 
in this study are part of the Mahantango Creek watershed. 
WE-38 Mahantango Creek is one of the experimental catch- 
ments of the Agricultural Research Service (ARS) of the 
U.S. Department of Agriculture (USDA). The total drainage 
area is 7.2 km 2 and it includes a 0.64-krrr subcatchment, 
WD-38, which has been the subject of intensive hydrologic 
research [e.g., Engman , 1974; Etigman and Rogowski , 
19 74 ]. 

The catchments experience a humid climate, typical for 
the northeastern United States, with approximately 1000 mm 
of annual precipitation, distributed uniformly throughout the 
year. The average rainfall during the month of July (the 
month of the machydro 90 experiment) is 104 mm [ Gburek , 
1977], The frequent rains keep the soils at or near field 
capacity, except near the surface. Even the fragipan soils 
near the stream, which have only limited moisture holding 
capacity and root penetration as a result of the confining 
layer, remain relatively moist because of this frequent rain- 
fall and their proximity to the high water table. Soil moisture 
is generally variable only to about the 100 cm depth [ Gburek , 
1977]. Below this depth soil moisture content remains nearly 
constant at near field capacity throughout the year. Based on 
measurements from more than 30 groundwater wells, Urban 
[1977] constructed a water table map for WE-38. The 
groundwater profile is highly correlated with topography. 
Because of the relatively high moisture content of the soils 
the actual evapotranspiration (ET) is at or near its potential 
rate much of the time. There are short periods, however, 
during which soil moisture deficit limits evapotranspiration. 
Gburek [1977] estimated that baseflow at WE-38 is about 
70% of yearly streamflow. 

3.2. Hydrologic Data 

The hydrologic data used in this study were gathered 
during the 12-day Multisensor Airborne Campaign 1990 
(machydro 90). The objective of this campaign is to study 
the role of soil moisture and its regulating effect on hydro- 
logic processes. The experiment was held from July 9 to July 
20, 1990. During this period the spatial and temporal distri- 
bution of precipitation was measured by 14 rain gages with 
0.2 mm accuracy at a sampling interval of 15 min. Figure 1 
shows the catchment average rainfall recorded during the 
12-day period. The observed accumulated rainfall is 66 mm. 
This rainfall fell during the first half of the experiment. The 
maximum observed rainfall intensity is about 20 mm/hour, a 
value low’er than the minimum saturated hydraulic conduc- 
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Fig. 1. (a) Areal average precipitation data for the 12-day 

period July 9-20, 1990 calculated from 15-min observations at 14 
rain gages, (b) Estimated potential evaporation based on 30-min 
measured net radiation and estimates based on weather reports (flat 
lines at the beginning and end of the graph). 

tivity reported for the soils in Mahantango Creek. Infiltration 
excess runoff is thus impossible during this period. 

Potential evaporation from July 1 1 to July 19 is calculated 
based on the Priestley-Taylor method using measured net 
radiation from a micrometeorological station. Radiation was 
measured at a sampling interval of 30 min. For July 9, 10, 
and 20, when the meteorological station was not operating, 
daily potential evaporation values of 5, 2.5, and 5 mm, 
respectively, are estimated based on weather reports. The 
total potential evaporation for the 12-day period is 45 mm. 
The temporal evolution of potential evaporation is shown in 
Figure 1. Based on these meteorological observations a 
strong drydown can be expected during the second half of 
the experiment. 

The temporal evolution of the phreatic surface was ob- 
served at four piezometers. The measured depth to the water 
table ranged from 1 to 4 m. Soil samples were taken along 
fixed transects, from which the gravimetric and volumetric 
moisture content in the top soil layer could be estimated. 
These data are net used in this study but will be analyzed and 
compared with the models’ output in future studies. Stream- 
flow is observed at the outlet of catchment WE-38 and 
subcatchment WD-38 (Figure 2). The level response at 180 
hours in catchment WE-38 is an artifact. No attempts were 
made to . correct for this anomalous observation since no 
detailed information about the data collection was available 
to the authors. During the 12-day experiment, it is estimated 
that only 1 or 2% of the total rainfall became direct runofi. 
Base flow tends to be low during summer months. During the 
experiment the volume of base flow was about equal to the 
direct runoff volume. The remaining rainfall observed during 
the experiment is either stored in the unsaturated and 
groundwater zones or is evaporated into the atmosphere. 
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Fig. 2. Observed streamflow data for the 12-day period July 9^20, 
1990 for (a) catchment WE-38 and ( b ) catchment WD-38. 


Figure 2 demonstrates the nonlinear response to rainfall of 
the two basins during the experiment. The initial storm 
response is almost negligible. The runoff production of the 
two basins during the last rainfall event is more pronounced 
and is dominated by the saturation excess mechanism. 

3.3. Model Parameters 

3.3.1. Soils. In WE-38 Mahantango Creek, 15 different 
soil types can be identified (Table 1). From Table 1 it can be 
seen that the spatial distribution of 0 S , 6 r , (i^ 5 ), and B(n) 

can be neglected for catchment WE-38, and therefore for 
subcatchment WD-38 as well. The following parameter 
values are used to characterize the soil water retention 
properties in WE-38 and WD-38: 9 S = 0.501, 8 r = 0.015, 
ip c = 0.21 m, ip s = 0.43 m, B = 0.211, and n = 1.29. The 
areal average^ value of saturated hydraulic conductivity at 
the surface K Sq is 0.062 m/hour. However, this parameter 
varies considerably from 0.036 to 0.090 m/hour and therefore 
is explicitly taken into account in both models. We refer to 
Rogowski et al . [1974] and Loague and Freeze [1985] for a 
detailed soil map of both WE-38 and WD-38. The value of 
specific storage S s used in (13) and (14) is 0.005/m. For the 
numerical model the parameter which controls the 

switching of evaporation boundary conditions, was set at a 
value low enough to ensure that actual evaporation remained 
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Fig. 3. Distribution of the local topographic index 10 In («/tan ft) 
for (a) catchment WE-38 and ( b ) catchment WD-38. 


at its potential rate throughout the simulation, in accordance 
with field observations. 

3.3.2. Topography . Catchment topography is repre- 
sented by a 30 x 30 m U.S. Geological Survey digital 
elevation model (DEM). From this DEM the topographic 
index In (tf/tan ft) can be determined for each grid square. 
Figure 3 shows the distribution of the topographic index in 
the two catchments. Sivapalan et al. [1987] used a three- 
parameter gamma distribution to model the topographic 
index: 


fiix) = 


1 

xU4>) 



exp 



(16) 


where x = In (a ! tan ft), <f> is a shape parameter, \ is a scale 
parameter, p is a shift parameter, and T( ) represents the 
gamma function. The parameter values for (16) applied to the 
two catchments were estimated by the method of moments 
and are given in Table 2. 

3.3.3. Numerical parameters . In the numerical model 
the vertical soil profile is discretized into six layers of 
thickness 5, 5, 12.5, 75, 127.5, and 275 cm, with the thinnest 
layers closest to the surface. This results in an impervious 
layer (the base of the catchment) at a depth of 5 m, running 
parallel to the land surface. The horizontal discretization is 
taken equal to the grid size in the digital elevation model, 
that is 30 m by 30 m. The spatial discretization of subcatch- 
ment WD-38 yields 3804 elements and 4935 nodes. The 


TABLE 2. Topographic Characteristics for Catchments WE-38 and WD-38 


Topographic Index, In {a ! tan ft) 


Coefficient of 

Catchment Mean Variance Variation 


Gamma Distribution Parameters 


<t> X M 


WE-38 4.03 2.88 0.42 5.651 0.714 0.0 

WD-38 4.37 2.48 0.36 T700 0.567 0.0 
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sensitivity of the simulated results to the grid resolution is an 
important issue, especially for areas in the catchment with 
large vertical pressure variations. Such areas would be near 
the streams and other runoff producing areas. Panic otii and 
Wood [1993] investigate the accuracy of the model to grid 
resolution and show that for grid sizes like those used in this 
study, good model simulations will be obtained. An adaptive 
time stepping algorithm is used in the numerical model, with 
a minimum time step size of 5 min and a maximum of 1 5 min 
(equal to the sampling interval of the rainfall data). During 
each time step the numerical solution is assumed to converge 
when the maximum change in nodal pressure head between 
nonlinear iterations is less than 6 cm. An average of four 
nonlinear iterations per time step is required for the simula- 
tions. 


3.4. Initial Conditions 

In the conceptual model the average depth to the water 
table z determines the initial catchment wetness conditions. 
For short-term simulations the output of the conceptual 
model is highly sensitive to this parameter, and it should 
therefore be carefully determined. Recently, Troch et al. 
[1993] developed a physically meaningful technique to esti- 
mate the effective depth to the water table, as a measure of 
the initial storage capacity of a basin. The analysis is based 
on Boussinesq’s hydraulic groundwater equation and uses 
streamflow measurements at the outlet of the basin. Bouss- 
inesq’s equation describes the water table height h(x , /) in 
the case of outflow into a stream channel from an idealized 
unconfined aquifer with width L placed on a horizontal 
impermeable layer: 


dh kd dh\ 

— = U— (17) 

dt n f dx y dx J 


where r represents the time since the start of the recession, 
jc is horizontal distance, h(x, t ) is the transient free ground- 
water surface profile, k is the hydraulic conductivity of the 
aquifer, and n € is the drainable or effective porosity (some- 
times referred to as specific yield). For small /, as the 
outflow rate q at x = 0 starts, the effect of the impermeable 
wall at x - L (representing the divide) is negligible and the 
solution to (17) can be taken to be the same as if L - «. 
Polubarinova-Kochina [1962] has presented an exact solu- 
tion for this case. The response of the aquifer water table to 
sudden drainage at t = 0 is not unlike a propagating wave. 
As soon as the wave reaches the end of the aquifer, at x = 
L, the small t solution is no longer valid. At this point, 
Boussinesq’s solution can be used. Boussinesq [1904] ob- 
tained an exact solution to the nonlinear differential equation 
by assuming that the initial water table has the form of an 
inverse incomplete beta function [see also Polubarinova- 
Kochina, 1962, pp. 515-517]. This solution is not valid for 
small t y when the aquifer is close to being fully saturated. 
Based on this large time solution, Troch et al. [1993] derived 
the following expression for the effective w'ater table height: 


h = 0.773D 


1 + 1.115 



' -i 

/ 


(18) 


where D represents depth to the bedrock. The catchment 
base flow' Q h can now be expressed as [ Troch et al., 1993tf] 


Q h = S.lllkiD - z) z D d L t (19) 


where k and D can be considered as catchment scale 
effective values of hydraulic conductivity and depth to the 
bedrock, respectively. D (i represents drainage density, and 
L, is the total length of the perennial channels. One can 
determine D. once k is know n, by defining a critical value of 
the base flow Q c . corresponding to a situation where it is 
assumed that the aquifers start to behave in accordance with 
Boussinesq’s solution: 

q =, 3.450 kD z D d L t (70) 

Equation ( 19) can then be used to estimate z from observed 
base flow at the start of the simulation period. 

To estimate the parameters in ( 19) and (20), a method for 
base flow analysis based on the following relationship is 
adopted: 


dQidt = <J>< Q) (21) 

Brutsaert and Nieber [1977] have shown that, for several 
solutions based on the Dupuit-Boussinesq hydraulic theory, 
d>( ) can be written in the form of a power function: 

dQidt = -a x Q h] (22) 


where a^ and b\ can be related to hydraulic and geomor- 
phologic characteristics of the basin. Using Boussinesq’s 
solution it follows immediately that [Brutsaert and Nieber , 
1977] 


O] 


4.804A J/: L, 
n e A ~' 2 


b j = 3/2 


(23) 


For the small time solution it can be shown that 


a i 


1.133 

kn f D y Lj 


(24) 


In Figure 4 we have plotted historical observed daily 
recession flow data for WE-38 and WD-38 versus its time 
derivative. Only uninterrupted recession flow data starting 
the second day after the cessation of rainfall are considered. 
For WE-38 (Figure 4 a) the data are fitted by a regression line 
(not shown on the figure) with slope 1.498 and a correlation 
coefficient of 0.88; the regression for WD-38 (Figure 4 b) has 
a slope of 1.35 and a correlation coefficient of 0.89. This 
suggests that the Dupuit-Boussinesq hydraulic theory holds: 
the observed slope is close to the theoretical slope of 1.5. 
Figure 4 shows a lower envelope with slope 1 .5 that excludes 
about 5 % of the data points. A lower envelope is chosen to 
eliminate other outflow components such as overland flow, 
interflow, channel drainage, and evaporation losses. How- 
ever, the exact position of this lower envelope is uncertain. 
Troch et al. [1993] suggest the use of a 5 or \0% lower 
envelope. In this study a lower envelope excluding 5% of the 
data points is adopted. The corresponding intercept values 
a f are 1.7861 x 10~ 6 (m 3 /s)’ l,: for WE-38 and 8.6792 x 
ICT 6 (m 3 /s) _1 ‘ for WD-38. 

The critical base flow value Q ( represents an upper limit 
for the applicability of the large time solution. Troch et a!. 
[19931 therefore suggest that this parameter should be esti- 
mated by the abscissa value of the intersection of the lower 
envelope curves with slopes 1.5 and 3 on a log-log diagram. 
They observed that, for a river basin in Belgium, the value of 
Q c is rather insensitive to the actual position ol these lines 
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Fig. 4. Log-log plot of — dQldt versus observed discharge Q for 
(a) catchment WE-38 and (b) catchment WD-38; also shown are the 
lower envelopes (solid lines) which exclude 5% of the data points. 


However, for the catchments under study, the slope of 3 is 
not apparent in the data. Therefore the critical value is 
estimated to be equal to the maximal observed base flow 
value in Figure 4 (0.500 m 3 /s and 0.050 m 3 /s, respectively). 

Table 3 summarizes the results of applying (23) and (20). 
Results are given for a range of n e values, namely, from 0.02 
to 0.07. These are reasonable values for the catchments 
under study [Freeze and Cherry , 1979, p. 61]. The geomor- 
phologic parameters D d and L t are estimated from the DEM 
data by means of an automated extraction algorithm [Band, 
1986]. This algorithm produces a mapping of stream chan- 
nels, ridges, and drainage basins. The total length of peren- 
nial channels, as defined by the blue lines on the topographic 
maps for WE-38 is about 12 km. which results in a drainage 
density of 1 .6 x 1 0 ~ 3 m “ 1 ; corresponding values for WD-38 
are 0.9 km and 1.5 x 10‘ 3 m" 1 . 

The observed base flow at the beginning of the experiment 
is 0.006 m 3 /s for WE-38 and 0.001 m 3 /s for WD-38. By 
means of (17) and using an average value of drainable 
porosity of 0.04, the effective depth to the water table z is of 
the order of 3.3 m for WE-38 and 2.3 m for WD-38. These 


TABLE 3. Estimated Aquifer Parameters for Catchments 
WE-38 and WD-38 


Drainable 

Porosity 

tQr 

WE-38 

= 0.500 m J /s> 

WD-38 
«? r = 0.050 

m-Vs) 

k , cm/h /), m 

k , cm/h 

D, m 

0.02 

53 

7.16 

134 

5.19 

0.03 

120 

4.78 

301 

3.46 

0.04 

213 

3.58 

536 

2.60 

0.05 

333 

2.87 

837 

2.08 

0.06 

480 

2.39 

1205 

1.73 

0.07 

653 

2.05 

1641 

1.49 


n 

lllWiL 
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Fig. 5. («) Distribution of the local water table depth at the start 
of the numerical simulation as calculated from (2). ( b ) Distribution 
of the local water table depth at the end of the numerical simulations 
for the 12-day period of the machydro 90 experiment. 




values are used as initial conditions in the conceptual water 
balance model. 

The initial conditions required for a transient simulation 
with the numerical model are nodal pressure head values. 
The initial heads are generated based on knowledge of the 
initial water table distribution given by (2) (see Figure 5a). 
The local water table depth is converted into a vertical 
pressure head distribution using the hydrostatic assumption. 
This means that if at a given location the depth to the water 
table is 3 m, the pressure head at the surface node is -3 m. 

3.5. Model Calibration 

Both the conceptual and the numerical model were cali- 
brated in order to preserve total runoff volume. During the 
simulation period the dominating runoff production mecha- 
nism is saturation excess. This is not surprising considering 
the rainfall characteristics during machydro 90 and the soil 
characteristics of the basin. The parameter / was used in 
both models as a fitting parameter. For different soil types 
and land use. Seven [1982] reports fitted values for the 
parameter / to observed soil hydraulic characteristics. The 
values range from 1.0 m" 1 to about 10. 0 m _l . For sandy 
loam to silt loam soils a typical value of 2.5 m ! is suggested. 
For one of the soil types in the catchment (Albrights; see 
Table 1). Rogowski et al. (1974] estimated this parameter to 
be of order 1 .2 m ~ l . In this study a catchment-wide effective 
value of / = 2.5 m~ l was found for the conceptual nodel, 
and a value of 1.05 m was used for the numerical model. 

It is noted that the optimized values of / correspond to 
observed values reported in the literature. 

The model parameters that w'ere estimated or inferred 
from data outside the 12 day simulation period are as 
follows: soil parameters described in section 3.3. 1 (estimated 
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from detailed soil survey data), topographic parameters 
described in section 3.3.2 (from LSGS DEM data) and initial 
water table depths described in section 3.4 (inferred from 
historical streamfiow' recession data). This leaves only J 
which was estimated during model calibration, as described 
above. 

4. Evaluation of thf Conceptual 
Water Balance Mo dll 

4.1. Spatial Distribution of Local Water Tabic Depth 

In the conceptual model the difference between the local 
water table depth and its areal average is a linear function of 
the combined topography-soil index In (aT c !T t tan /3), as 
expressed by (2). Based on field evidence. Wood et al. [ 1990] 
concluded that the deviation of the topographic variable 
from its expected value A is far greater than the standard 
deviation of the local values of the transmissivity coefficient. 
Therefore variability in the transmissivity coefficient will 
have a relatively small effect on the predicted patterns of 
water table depths in the conceptual model compared to the 
effect of topographic variability. The observations of Urban 
[1977] support this result for Mahantango Creek. Figure 5a 
shows the initial distribution of local water table depth for 
the WD-38 catchment calculated based on (2). This distribu- 
tion is used in the initial conditions for both the conceptual 
and the numerical model. During simulations with the nu- 
merical model the water table profile, as a function of space 
( jc and y coordinates) and time, is calculated based on the 
three-dimensional equations governing flow in porous me- 
dia. At the end of the 12-day simulation period which 
involved hydrologic fluxes of significant magnitude (66 mm 
of rainfall and 45 mm of potential evaporation), the numer- 
ical model preserves, to a certain extent, the initial distribu- 
tion of the water table depth. This suggests that the estimate 
of the initial conditions based on (2) is consistent with 
groundwater hydraulics and that the use of (2) in a basin 
scale water balance model is reasonable. 

To test the effect of initial conditions, longer numerical 
simulations for the WD-38 catchment were performed. Dur- 
ing the first 100 days of a 400-day simulation, we apply 
zero-flux boundary conditions at the surface (no precipita- 
tion and no evaporation), in order to allow redistribution 
within the saturated zone. For the next 150 days a constant 
rainfall rate of 0.1 mm/hour is applied as the surface bound- 
ary condition. This is followed by a second zero-flux bound- 
ary condition period of 150 days. The distribution of the 
water table depth at the end of the first 100 days is shown in 
Figure 6a, while the distribution at the end of the simulation 
period for the same catchment is given in Figure 6b. After 
the first 100-day period the distribution of the local water 
table depth closely resembles the initial distribution. This 
means that, even after a longer simulation period, the model 
shows no tendency to drift away from the initial conditions 
set by (2). At the end of the 400-day simulation, however, the 
distribution tends to become more uniform. We also observe 
a shift in the areal average water table depth due to the large 
amount of rainfall (360 mm) during this period. We conclude 
from these tests that the linear relationship between water 
table depth and topographic index given by (2) appears to be 
a reasonable assumption except when the catchment is 
highly stressed or far removed from its steady or equilibrium 
state as may occur during prolonged periods of rainfall. 
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Fig. 6. (a) Distribution of the local water table depth after 100 

days of a 400-day numerical simulation, with zero-flux boundary 
conditions at the surface for the first 100 days. ( b) Distribution of the 
local water table depth at the end of the 400-day numencal simula- 
tion, with total rainfall of 360 mm. 

4.2. Temporal Evolution of Water Table Depth 

During a rainstorm the conceptual model does not update 
the depth to the water table. Immediately after cessation of 
rainfall a new value of mean water table depth is calculated 
taking into account the infiltrated and drained volume from 
the previous storm. This aspect of the model is apparent in 
Figure la. Line 1 in Figure la shows the variation of areal 
average water table depth about its mean value, calculated 
with the conceptual model for WD-38. The mean water table 
depth during the simulation is 2.13 m. During the dry down 
at the end of the experiment the evolution of mean water 
table depth is controlled by evaporative and drainage losses. 
Line 2 in Figure la shows the variation of areal average 
water table depth about its mean value, as computed with 
the numerical model for WD-38. The calculated mean value 
is 2.85 m. The updating of the water table is now controlled 
by the percolation and a much smoother curve is obtained. 

The variation of local water table depth with respect to its 
local mean as computed by the conceptual model is similar 
to line 1 in Figure la. This is due to the fact that the soil 
characteristics which are controlling the updating of the 
water table are assumed to be spatially invariant and that 
runoff production is dominated by the saturation excess 
mechanism. Therefore we can compare the areal average 
simulation results for the conceptual model to observed 
variations in water table for the catchment. Line 3 in Figure 
lb shows the simulated variation in water table depth by- 
means of the numerical mode! for a site along the transect 
shown in Figure 8. The mean water table depth for this site 
during the experiment is 1.43 m. Comparing line 3 in Figure 
lb with line 2 in Figure 7a we can see that the variation of 
areal average water table depth calculated with the numeri- 
cal model reflects the variation for sites in the catchment 
exhibiting deeper water table depths. Therefore it is inter- 
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Fig. 7. (a) Line 1: variation about the mean value (2.13 m) of 

the areal average water table depth (WTD) calculated with the 
conceptual model for the 12-day period; line 2: variation about the 
mean value (2.85 m) of the areal average WTD calculated with the 
numerical model for the 12-day period; line 3: observed variation of 
WTD in a deep well (mean value 3.63 m); and line 4: observed 
variation of WTD in a shallow well (mean value 0.91 m). This line is 
shifted from the zero mean to overlap with line 1 . (b) Numerically 
simulated evolution of local water table depth at four locations along 
a transect; line 1: mean WTD value 0.32 m (node closest to the 
channel); line 2: mean WTD value 1.13 m; line 3: mean WTD value 
1.43 m; and line 4; mean WTD value 1.92 m (node farthest from the 
channel). 


esting to compare the areal average simulation results for the 
numerical model to observed variations in deeper piezome- 
ters. Line 4 in Figure la shows the observed variation in 
depth to water table for a shallow piezometer (mean value 
during the experiment: 0.91 m) and line 3 shows the ob- 
served variation for a deeper well (mean value: 3.63 m). The 
irregular fluctuation of the shallow' well observations is 



Fig. 8. Elevation image of subcatchment WD-38 (30 x 30 m grid 
resolution) showing stream (heavily shaded pixels) and location of 
the transect of four surface nodes selected for vertical profile output 
(unshaded pixels). 


probably due to evapotranspiration losses. The groundwater 

vdrograph. as observed by the shallow piezometer, shows 
more or less the same temporal evolution as simulated by the 
conceptual model. During the last rainfall period the water 
balance model predicts a rise in w ater table of about 2-3 cm 
and the observations show a comparable evolution. After the 
rainfall the simulated drop in water table is in good agree- 
ment with the response in the shallow well. After a few days, 
however, the conceptual model overestimates the rate of 
decline of the water table. The evolution of mean w-ater table 
depth as calculated by the numerical model and the obser- 
vation for the same period in the deeper well are in remark- 
ably good agreement. The absolute mean values arc different 
(2.85 and 3.63 m, respectively) but the general trend in lines 
2 and 3 is very close. The diurnal variation in the observation 
in the deep well during the last days of the experiment 
(which is also exhibited by other piezometers, not shown in 
Figure 7) is explained by diurnal variations in atmospheric 
pressure. Similar results were obtained using data observed 
in the other piezometers. 

These results seem to indicate that the simplifying as- 
sumption in the conceptual model concerning the temporal 
evolution of the water table holds for those parts of the 
catchment with a shallow water table. This could be ex- 
plained by the fact that the time delay between infiltration 
and percolation to the water table is small for these sites. To 
further test this hypothesis we computed, from the numeri- 
cal simulation of WD-38, the temporal evolution of the water 
table for different locations along a hillslope transect. The 
four unshaded pixels in Figure 8 show the location of the 
surface nodes selected for detailed vertical profile output. 
The results for the 12-day simulation period are given in 
Figure lb. It is clear from this graph that infiltration and 
percolation effects are more damped uphill (line 4), where 
the water table is deeper. The range of variation of depth to 
the water table close to the channel (line 1) compares 
reasonably well to the one predicted by the conceptual 
model. It has a mean value of 0.30 m. 

It is interesting to refer to earlier research concerning 
groundwater level fluctuations performed in the same basin. 
Urban [1977] reports groundwater recharge events for three 
wells in September 1973. For the shallow well (water table at 
about 1 m below the land surface) a similar response to 
rainfall, as observed in the shallow well during machydro 90, 
is reported. For a comparable rainfall event the variation is 
a few centimeters. However, the response measured in the 
deep well (water level at about 18 m below the land surface) 
is much more pronounced than that observed during machy- 
dro 90, with a variation of the order of 1 m. The measured 
water table in this well is in the outcrop of a sandstone. The 
sandstone is covered with coarse gravelly colluvial soils. 
These geologic conditions are not represented in the numer- 
ical simulations reported here. 

4.3. Soil Moist it re Profi I e s 

The conceptual model calculates unsaturated zone storage 
capacity based on the hydraulic equilibrium assumption. 
This assumption also afl'ects the local infiltration parameters 
needed in Philip's equation. Figure 9 shows simulated pres- 
sure head profiles produced by the numerical model for 
catchment WD-38. The locations of these vertical profiles in 
the catchment are those shown in Figure 8. We can see that 
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Fig. 9. Vertical pressure head profiles for five different time steps during the 12-dav numerical simulation and at the 
four locations along the transect shown in Figure 8: (a) node farthest from the channel; ( d ) node closest to the channel, 


in the upper 1 m of the unsaturated zone the pressure head 
profiles are not consistent with the hydraulic equilibrium 
assumption. This fact is detected for the pixels furthest from 
the stream where the water table is deepest. Below 1 m 
depth the simulated profiles do not deviate significantly from 
hydraulic equilibrium. Based on field evidence collected at 
the Mahantango Creek watershed, Gburek [1977] concluded 
that variation in soil moisture is limited to the upper 1 m 
layer, and that below this depth soil moisture content 
remains nearly constant and near field capacity throughout 
the year. 

These results suggest that we can improve the conceptual 
model by using two layers instead of one to model the 
unsaturated zone. The upper layer in such a two-layer model 
can also incorporate root zone processes. This extension has 
been made to the conceptual model and is being tested on the 
First 1SLSCP Field experiment (FIFE) data set from Kansas 
[ Famiglietti , 1992; J, S. Famiglietti and E. F. Wood, Aggre- 
gation and scaling of spatially variable hydrological pro- 
cesses, 2, A catchment scale model of water an energy 
balance, submitted to Water Resources Research , 1993]. 

4.4. Base Flow Recession Characteristics 

The parameters a \ and b ] of (22) in the case of an 
exponential subsurface saturated soil water store are given 
by 

a \ = f IA b x = 2 (25) 

In contrast to observations for catchment WE-38 and sub- 
catchment WD-38 (Figure 4) and to Boussinesq's hydraulic 
groundwater theory (equations (23)), the predicted subsur- 
face flow contributions from the conceptual model will yield 
a slope of 2 on a log KdQldt) versus log Q diagram (see line 


2 in Figure 10). We believe that this is not an unreasonable 
slope value for characterizing base flow recession in steeper 
catchments, where the influence of the hillslopes on ground- 
water flow is significant, at least during the initial stage of a 
recession period, and thus where Boussinesq's theory is not 
valid [Zecharias and Brutsaert , 1988]. For the catchments 
used in this study, however, topographic effects on base flow 
recession characteristics are negligible and base flow obser- 
vations correspond to Boussinesq's groundwater equation 
and therefore give a slope close to 1.5 on the log (dQ/dt) 
versus log Q diagram. For catchments with mild slopes the 
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Fig. 10. Theoretical versus .simulated base flow recession char- 
acteristics. X: log-log plot of -JQ dt versus discharge Q for 
numerically simulated base flow contribution; line I: least squares fit 
to the data points (slope 1.57); and line 2; theoretical base flow 
characteristic from the conceptual model (slope 2). 
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base flow equation in the conceptual model should be 
modified in order to fully represent the recession dynamics. 

Since the numerical model solves the three-dimensional 
Richards equation in an isotropic soil matrix it is not 
surprising that the calculated subsurface flow contributions 
(plotted points in Figure 10) during the drydown period in the 
simulation run for subcatchment WD-38 behave as predicted 
by Boussinesq’s theory. In fact, the results from the numer- 
ical model yield a slope of 1.57 in Figure 10 (line 1). 

Equations (25) offer the possibility of estimating the model 
parameter/ based on recession observations. By imposing a 
slope of 2 on the data displayed in a log ( dQ/dt ) versus log 
Q diagram and by taking a lower envelope excluding 5% of 
the data, the value of the intercept can be calculated. 
Through the first of equations (25) a corresponding value of 
/ can then be estimated. In (25) the drainage area plays the 
role of a scaling factor. A similar technique can be used to 
determine the parameter/ in the numerical model. In future 
studies, different methods of calibrating both the numerical 
and conceptual models will be compared. 

5. Summary 

The evaluation of a conceptual catchment scale water 
balance model is presented. This evaluation is based on a 
comparison of simulation results from the conceptual model 
with results from a three-dimensional physically based nu- 
merical model and with field observations. The conceptual 
model relies on a topographic index to predict saturation 
excess runoff and on Philip's equation to predict infiltration 
excess runoff. The numerical model is based on the three- 
dimensional transient Richards equation. The study is car- 
ried out at two basin scales: a 7,2-km 2 catchment and a 
0.64-km : subcatchment. These catchments are located in the 
North Appalachian ridge and valley region of eastern Penn- 
sylvania. Detailed hydrologic data were collected for both 
catchments during the 12-day machydro 90 experiment. 

Simplifying assumptions in the conceptual model concern- 
ing the spatial distribution and temporal evolution of water 
table depth, the distribution of soil moisture and pressure 
head in the unsaturated zone, and the characteristics of base 
flow recession are discussed. The hypothesis about the 
spatial and temporal evolution of local water table depth is 
tested by comparing the distribution of the topographic 
index with the distributions of the water table depth gener- 
ated by the numerical model at the end of a 12-day simula- 
tion period. It is found that the use of a linear relationship 
between the local water table depth and the topographic 
index is reasonable in a basin scale water balance model. To 
test the effect of initial conditions, longer test runs with the 
numerical model are performed. Piezometric observations 
are compared to simulated groundwater dynamics. Close 
agreement between simulation results and field observations 
is obtained for shallow well observations. The hydraulic 
equilibrium assumption for the vertical distribution of soil 
moisture in the conceptual model is tested by examining the 
moisture profiles obtained from the numerical model. The 
results suggest that the conceptual model can be improved 
by using two layers to model the unsaturated zone. The 
characteristics of base flow recession for the conceptual and 
numerical models are compared with analytical solutions to 
Boussinesq's hydraulic equation and with observations. 
Based on this comparison we suggest that strearnfiow reces- 


sion data can be used to calibrate both the conceptual model 
and the numerical model. 

Further tests are required to confirm some of the findings 
discussed in this paper. Based on the 12-day simulation 
period for subcatchment WD-38 we conclude that the distri- 
bution of the topographic index is a reasonable measure for 
the local water table depth in Mahantango Creek. Detailed 
hydrologic data, as used in this study but for a longer period, 
are necessary to further justify this conclusion. The effect of 
catchment topography on the base flow recession character- 
istics, as suggested in section 4.4, should be further investi- 
gated. Other simplifying assumptions in the conceptual 
model, such as the calculation of soil moisture content, 
should be tested by comparing the simulation results with 
field observations and with remotely sensed information. A 
study comparing simulated soil moisture maps with soil 
moisture maps generated from remotely sensed information 
for WE-38 and WD-38 Mahantango Creek and based on 
machydro 90 data is in progress. This work will also involve 
extending some of the numerical simulation tests to WE-38 
and to other catchments. 
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